﻿using System;
using System.Collections.Generic;
using System.Linq;
using System.Web;

/// <summary>
/// GIS 全坐标转换类
/// 包括  WGS1984
///       火星坐标
///       百度
///       Web Mercator
///       距离计算
/// </summary>
/// <remarks>
/// http://www.soshow.org/html/25/n-52025.html
/// http://www.csdn.cc/post/2014-01-20/40060829205
/// http://www.oschina.net/code/snippet_260395_39205
///由php坐标转换改变 
/// </remarks>
public class MapGisCoodConverter
{

    private const double PI = 3.14159265358979324;
    private const double x_pi = 3.14159265358979324 * 3000.0 / 180.0;

    /// <summary>
    /// WGS-84 to GCJ-02
    /// </summary>
    /// <param name="wgsLat"></param>
    /// <param name="wgsLon"></param>
    /// <returns></returns>
    public static double[] gcj_encrypt(double wgsLat, double wgsLon)
    {
        if (outOfChina(wgsLat, wgsLon))
            return new double[] { wgsLat, wgsLon };
        double[] d = delta(wgsLat, wgsLon);
        return new double[] { wgsLat + d[0], wgsLon + d[1] };
    }
    /// <summary>
    /// GCJ-02 to WGS-84
    /// </summary>
    /// <param name="gcjLat"></param>
    /// <param name="gcjLon"></param>
    /// <returns></returns>
    public static double[] gcj_decrypt(double gcjLat, double gcjLon)
    {
        if (outOfChina(gcjLat, gcjLon))
            return new double[] { gcjLat, gcjLon };

        double[] d = delta(gcjLat, gcjLon);
        return new double[] { gcjLat - d[0], gcjLon - d[1] };
    }
    /// <summary>
    /// GCJ-02 to WGS-84 exactly
    /// </summary>
    /// <param name="gcjLat"></param>
    /// <param name="gcjLon"></param>
    /// <returns></returns>
    public static double[] gcj_decrypt_exact(double gcjLat, double gcjLon)
    {
        double initDelta = 0.01;
        double threshold = 0.000000001;
        double dLat = initDelta;
        double dLon = initDelta;
        double mLat = gcjLat - dLat;
        double mLon = gcjLon - dLon;
        double pLat = gcjLat + dLat;
        double pLon = gcjLon + dLon;
        double wgsLat = 0;
        double wgsLon = 0;
        int i = 0;
        while (true)
        {
            wgsLat = (mLat + pLat) / 2;
            wgsLon = (mLon + pLon) / 2;
            double[] tmp = gcj_encrypt(wgsLat, wgsLon);
            dLat = tmp[0] - gcjLat;
            dLon = tmp[1] - gcjLon;
            if ((Math.Abs(dLat) < threshold) && (Math.Abs(dLon) < threshold))
                break;

            if (dLat > 0) pLat = wgsLat; else mLat = wgsLat;
            if (dLon > 0) pLon = wgsLon; else mLon = wgsLon;

            if (++i > 10000) break;
        }
        //console.log(i);
        return new double[] { wgsLat, wgsLon };
    }
    /// <summary>
    /// GCJ-02 to BD-09
    /// </summary>
    /// <param name="gcjLat"></param>
    /// <param name="gcjLon"></param>
    /// <returns></returns>
    public static double[] bd_encrypt(double gcjLat, double gcjLon)
    {
        double x = gcjLon;
        double y = gcjLat;
        double z = Math.Sqrt(x * x + y * y) + 0.00002 * Math.Sin(y * x_pi);
        double theta = Math.Atan2(y, x) + 0.000003 * Math.Cos(x * x_pi);
        double bdLon = z * Math.Cos(theta) + 0.0065;
        double bdLat = z * Math.Sin(theta) + 0.006;
        return new double[] { bdLat, bdLon };
    }
    /// <summary>
    /// BD-09 to GCJ-02
    /// </summary>
    /// <param name="bdLat"></param>
    /// <param name="bdLon"></param>
    /// <returns></returns>
    public static double[] bd_decrypt(double bdLat, double bdLon)
    {
        double x = bdLon - 0.0065;
        double y = bdLat - 0.006;
        double z = Math.Sqrt(x * x + y * y) - 0.00002 * Math.Sin(y * x_pi);
        double theta = Math.Atan2(y, x) - 0.000003 * Math.Cos(x * x_pi);
        double gcjLon = z * Math.Cos(theta);
        double gcjLat = z * Math.Sin(theta);
        return new double[] { gcjLat, gcjLon };
    }
    /// <summary>
    ///WGS-84 to Web mercator
    ///mercatorLat -> y mercatorLon -> x
    ///
    /// </summary>
    /// <param name="wgsLat"></param>
    /// <param name="wgsLon"></param>
    /// <returns></returns>
    public static double[] mercator_encrypt(double wgsLat, double wgsLon)
    {
        double x = wgsLon * 20037508.34 / 180;
        double y = Math.Log(Math.Tan((90 + wgsLat) * PI / 360)) / (PI / 180);
        y = y * 20037508.34 / 180;
        return new double[] { y, x };
        /*
        if ((abs(wgsLon) > 180 || abs(wgsLat) > 90))
            return NULL;
        x = 6378137.0 * wgsLon * 0.017453292519943295;
        a = wgsLat * 0.017453292519943295;
        y = 3189068.5 * log((1.0 + sin(a)) / (1.0 - sin(a)));
        return  new double[]{y,  x);
        //*/
    }
    /// <summary>
    /// Web mercator to WGS-84
    /// mercatorLat -> y mercatorLon -> x
    /// </summary>
    /// <param name="mercatorLat"></param>
    /// <param name="mercatorLon"></param>
    /// <returns></returns>
    public static double[] mercator_decrypt(double mercatorLat, double mercatorLon)
    {
        double x = mercatorLon / 20037508.34 * 180;
        double y = mercatorLat / 20037508.34 * 180;
        y = 180 / PI * (2 * Math.Atan(Math.Exp(y * PI / 180)) - PI / 2);
        return new double[] { y, x };
        /*
        if (abs(mercatorLon) < 180 && abs(mercatorLat) < 90)
            return NULL;
        if ((abs(mercatorLon) > 20037508.3427892) || (abs(mercatorLat) > 20037508.3427892))
            return NULL;
        a = mercatorLon / 6378137.0 * 57.295779513082323;
        x = a - (floor(((a + 180.0) / 360.0)) * 360.0);
        y = (1.5707963267948966 - (2.0 * atan(exp((-1.0 * mercatorLat) / 6378137.0)))) * 57.295779513082323;
        return  new double[]{y,  x};
        //*/
    }

    /// <summary>
    /// 计算两点间的距离 单位(米)
    ///  two point's distance
    /// </summary>
    /// <param name="latA"></param>
    /// <param name="lonA"></param>
    /// <param name="latB"></param>
    /// <param name="lonB"></param>
    /// <returns> 单位(米)</returns>
    public static double distance(double latA, double lonA, double latB, double lonB)
    {
        double earthR = 6371000;
        double x = Math.Cos(latA * PI / 180) * Math.Cos(latB * PI / 180) * Math.Cos((lonA - lonB) * PI / 180);
        double y = Math.Sin(latA * PI / 180) * Math.Sin(latB * PI / 180);
        double s = x + y;
        if (s > 1) { s = 1; }
        if (s < -1) { s = -1; }
        double alpha = Math.Acos(s);
        double distance = alpha * earthR;
        return distance;
    }

    /// <summary>
    /// 
    /// </summary>
    /// <param name="lat"></param>
    /// <param name="lon"></param>
    /// <returns></returns>
    private static double[] delta(double lat, double lon)
    {
        // Krasovsky 1940
        //
        // a = 6378245.0, 1/f = 298.3
        // b = a * (1 - f)
        // ee = (a^2 - b^2) / a^2;
        double a = 6378245.0;//  a: 卫星椭球坐标投影到平面地图坐标系的投影因子。
        double ee = 0.00669342162296594323;//  ee: 椭球的偏心率。
        double dLat = transformLat(lon - 105.0, lat - 35.0);
        double dLon = transformLon(lon - 105.0, lat - 35.0);
        double radLat = lat / 180.0 * PI;
        double magic = Math.Sin(radLat);
        magic = 1 - ee * magic * magic;
        double sqrtMagic = Math.Sqrt(magic);
        dLat = (dLat * 180.0) / ((a * (1 - ee)) / (magic * sqrtMagic) * PI);
        dLon = (dLon * 180.0) / (a / sqrtMagic * Math.Cos(radLat) * PI);
        return new double[] { dLat, dLon };
    }
    /// <summary>
    /// 判断是否中国范围
    /// </summary>
    /// <param name="lat">纬度</param>
    /// <param name="lon">经度</param>
    /// <returns>返回true,false</returns>
    public static bool outOfChina(double lat, double lon)
    {
        if (lon < 72.004 || lon > 137.8347)
            return true;
        if (lat < 0.8293 || lat > 55.8271)
            return true;
        return false;
    }

    /// <summary>
    /// 转换纬度
    /// </summary>
    /// <param name="x"></param>
    /// <param name="y"></param>
    /// <returns></returns>
    private static double transformLat(double x, double y)
    {
        double ret = -100.0 + 2.0 * x + 3.0 * y + 0.2 * y * y + 0.1 * x * y + 0.2 * Math.Sqrt(Math.Abs(x));
        ret += (20.0 * Math.Sin(6.0 * x * PI) + 20.0 * Math.Sin(2.0 * x * PI)) * 2.0 / 3.0;
        ret += (20.0 * Math.Sin(y * PI) + 40.0 * Math.Sin(y / 3.0 * PI)) * 2.0 / 3.0;
        ret += (160.0 * Math.Sin(y / 12.0 * PI) + 320 * Math.Sin(y * PI / 30.0)) * 2.0 / 3.0;
        return ret;
    }
    /// <summary>
    /// 转换经度
    /// </summary>
    /// <param name="x"></param>
    /// <param name="y"></param>
    /// <returns></returns>
    private static double transformLon(double x, double y)
    {
        double ret = 300.0 + x + 2.0 * y + 0.1 * x * x + 0.1 * x * y + 0.1 * Math.Sqrt(Math.Abs(x));
        ret += (20.0 * Math.Sin(6.0 * x * PI) + 20.0 * Math.Sin(2.0 * x * PI)) * 2.0 / 3.0;
        ret += (20.0 * Math.Sin(x * PI) + 40.0 * Math.Sin(x / 3.0 * PI)) * 2.0 / 3.0;
        ret += (150.0 * Math.Sin(x / 12.0 * PI) + 300.0 * Math.Sin(x / 30.0 * PI)) * 2.0 / 3.0;
        return ret;
    }

}

/*
 * 
http://www.soshow.org/html/25/n-52025.html
 解决GPS偏移代码
发布: 2013-11-13 09:46 |  作者: 閑闪互动 |   |  查看: 71次 

WGS-84 坐标系，是 GPS 系统所采用的坐标系。一切正常工作的 GPS 或 GPS 芯片所返回的坐标值都是这个坐标系下的数值。Google 地图采用的卫星图也是按照这个坐标系摆放的。

GCJ-02 坐标系，是我天朝政府搞出来的加密坐标系，也常常被称为 “火星坐标系”。包括（但可能不限于）高德地图在内的国内地图服务商采用它来绘制地图。Apple、Google 等国外公司在其道路地图中使用的也是高德的数据。BD-09 坐标系则是百度地图专用的坐标系。

在开发一些 LBS 应用时，如果不加处理，很容易出现几种形式的地图之间出现偏移的情况。因此在这几个坐标系之间进行转换非常重要。以下代码就是网络上泄露出的从 WGS-84 转 GCJ-02 的算法，以及 GCJ-02 与 BD-09 的互相转换算法。

另外，代码中的从 GCJ-02 到 WGS-84 的逆变换是我用迭代法实现的。实验证明，由于原变换的局部线性性非常好，只要迭代两次（调用原变换三次），经纬度就可以达到大约 1e-8 度的精度（相当于毫米级精度）。所以运行速度不会太慢，对于一般应用足矣。

全部代码也可以在 GitHub 下载：https://github.com/fengzee-me/ChinaMapShift。代码和其编译出的 libChinaMapShift.a 库文件直接适用于 iOS 环境下，需要在其他平台下运行时，可以自行编译（源代码完全是原生 C 代码，只需要 math.h 库）。


[cpp]view plaincopy
//  
//  ChinaMapShift.h  
//  ChinaMapShift  
//  
//  Most code created by someone anonymous.  
//  transformFromGCJToWGS() added by Fengzee (fengzee@fengzee.com).  
//  
   
#ifndef ChinaMapShift_ChinaMapShift_h  
#define ChinaMapShift_ChinaMapShift_h  
   
typedef struct {  
    double lng;  
    double lat;  
} Location;  
   
Location transformFromWGSToGCJ(Location wgLoc);  
Location transformFromGCJToWGS(Location gcLoc);  
Location bd_encrypt(Location gcLoc);  
Location bd_decrypt(Location bdLoc);  
   
#endif  


[cpp]view plaincopy
//  
//  ChinaMapShift.mm  
//  ChinaMapShift  
//  
//  Most code created by someone anonymous.  
//  transformFromGCJToWGS() added by Fengzee (fengzee@fengzee.com).  
//  
  
#include "ChinaMapShift.h"  
#include <math.h>  
  
inline Location  
LocationMake(double lng, double lat) {  
    Location loc; loc.lng = lng, loc.lat = lat; return loc;  
}  
  
///  
///  Transform WGS-84 to GCJ-02 (Chinese encrypted coordination system)  
///  
  
const double pi = 3.14159265358979324;  
  
//  
// Krasovsky 1940  
//  
// a = 6378245.0, 1/f = 298.3  
// b = a * (1 - f)  
// ee = (a^2 - b^2) / a^2;  
const double a = 6378245.0;  
const double ee = 0.00669342162296594323;  
  
int outOfChina(double lat, double lon)  
{  
    if (lon < 72.004 || lon > 137.8347)  
        return 1;  
    if (lat < 0.8293 || lat > 55.8271)  
        return 1;  
    return 0;  
}  
  
double transformLat(double x, double y)  
{  
    double ret = -100.0 + 2.0 * x + 3.0 * y + 0.2 * y * y + 0.1 * x * y + 0.2 * sqrt(x > 0 ? x:-x);  
    ret += (20.0 * sin(6.0 * x * pi) + 20.0 *sin(2.0 * x * pi)) * 2.0 / 3.0;  
    ret += (20.0 * sin(y * pi) + 40.0 * sin(y / 3.0 * pi)) * 2.0 / 3.0;  
    ret += (160.0 * sin(y / 12.0 * pi) + 320 * sin(y * pi / 30.0)) * 2.0 / 3.0;  
    return ret;  
}  
  
double transformLon(double x, double y)  
{  
    double ret = 300.0 + x + 2.0 * y + 0.1 * x * x + 0.1 * x * y + 0.1 * sqrt(x > 0 ? x:-x);  
    ret += (20.0 * sin(6.0 * x * pi) + 20.0 * sin(2.0 * x * pi)) * 2.0 / 3.0;  
    ret += (20.0 * sin(x * pi) + 40.0 * sin(x / 3.0 * pi)) * 2.0 / 3.0;  
    ret += (150.0 * sin(x / 12.0 * pi) + 300.0 * sin(x / 30.0 * pi)) * 2.0 / 3.0;  
    return ret;  
}  
  
Location transformFromWGSToGCJ(Location wgLoc)  
{  
    Location mgLoc;  
    if (outOfChina(wgLoc.lat, wgLoc.lng))  
    {  
        mgLoc = wgLoc;  
        return mgLoc;  
    }  
    double dLat = transformLat(wgLoc.lng - 105.0, wgLoc.lat - 35.0);  
    double dLon = transformLon(wgLoc.lng - 105.0, wgLoc.lat - 35.0);  
    double radLat = wgLoc.lat / 180.0 * pi;  
    double magic = sin(radLat);  
    magic = 1 - ee * magic * magic;  
    double sqrtMagic = sqrt(magic);  
    dLat = (dLat * 180.0) / ((a * (1 - ee)) / (magic * sqrtMagic) * pi);  
    dLon = (dLon * 180.0) / (a / sqrtMagic * cos(radLat) * pi);  
    mgLoc.lat = wgLoc.lat + dLat;  
    mgLoc.lng = wgLoc.lng + dLon;  
  
    return mgLoc;  
}  
  
///  
///  Transform GCJ-02 to WGS-84  
///  Reverse of transformFromWGSToGC() by iteration.  
///  
///  Created by Fengzee (fengzee@fengzee.com) on 13-11-9.  
///  
Location transformFromGCJToWGS(Location gcLoc)  
{  
    Location wgLoc = gcLoc;  
    Location currGcLoc, dLoc;  
    while (1) {  
        currGcLoc = transformFromWGSToGCJ(wgLoc);  
        dLoc.lat = gcLoc.lat - currGcLoc.lat;  
        dLoc.lng = gcLoc.lng - currGcLoc.lng;  
        if (fabs(dLoc.lat) < 1e-7 && fabs(dLoc.lng) < 1e-7) {  // 1e-7 ~ centimeter level accuracy  
          // Result of experiment:  
          //   Most of the time 2 iterations would be enough for an 1e-8 accuracy (milimeter level).  
          //  
            return wgLoc;  
        }  
        wgLoc.lat += dLoc.lat;  
        wgLoc.lng += dLoc.lng;  
    }  
  
    return wgLoc;  
}  
  
///  
///  Transform GCJ-02 to BD-09  
///  
const double x_pi = 3.14159265358979324 * 3000.0 / 180.0;  
Location bd_encrypt(Location gcLoc)  
{  
    double x = gcLoc.lng, y = gcLoc.lat;  
    double z = sqrt(x * x + y * y) + 0.00002 * sin(y * x_pi);  
    double theta = atan2(y, x) + 0.000003 * cos(x * x_pi);  
    return LocationMake(z * cos(theta) + 0.0065, z * sin(theta) + 0.006);  
}  
  
///  
///  Transform BD-09 to GCJ-02  
///  
Location bd_decrypt(Location bdLoc)  
{  
    double x = bdLoc.lng - 0.0065, y = bdLoc.lat - 0.006;  
    double z = sqrt(x * x + y * y) - 0.00002 * sin(y * x_pi);  
    double theta = atan2(y, x) - 0.000003 * cos(x * x_pi);  
    return LocationMake(z * cos(theta), z * sin(theta));  
} 


*/